Dynamical Breaking of Charge Neutrality in Intrinsic Josephson Junctions: Common 
Origin for Microwave Resonant Absorptions and Multiple-branch Structures in the 

I-V Characterisitcs 
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We demonstrate that both microwave resonant absorptions and multiple-branch structures in the 
I-V characterisitcs observed in intrinsic Josephson junctions (IJJ's) are caused by dynamical breaking 
of charge neutrality (DBCN) inside the atomic-scale superconducting layer. The Lagranginan for 
the time-dependent Lawrence-Doniach model incorporating the effect of the DBCN is proposed, and 
the longitudinal collective Josephson plasma mode is proved to exist based on the Lagrangian. On 
the other hands, the branching behaviors in the I-V curves are almost completely reproduced by 
careful numerical simulations for the model equation derived from the Lagragian. 

PACS numbers: 74.25.Fy, 74.50.-Fr, 74.80.Dm 



Intrinsic Josephson effects (IJE's) in highly anisotropic 
high temperature superconductors (HTSC's) have at- 
tracted growing interests , , and recent extensive ex- 
perimental studies have discovered two remarkable phe- 
nomena for IJE's, i.e., the microwave resonant absorp- 
tions (MRA's) H] and the multiple-branch structures 
(MBS's) in I-V characterisitics along the c-axis Q|. 

In Bi-2212 single crystals(SC's), the MRA has been 
observed for both longitudinal and transverse microwave 
configurations in cavities, and the observed dispersion re- 
lations have clearly been distinguished from one another 
for those configurations ||]. The MRA in the longitudi- 
nal configuration is caused by the excitaion of the lon- 
gitudinal Josephson plasma mode propagating along the 
c-axis, and the existence of such the charge fluctuation 
mode indicates that the charge neutrality can be dynam- 
ically broken inside the superconducting layer under the 
application of an AC electric field along the c-axis. Thus, 
a novel coupling between junctions due to the DBCN has 
been suggested and new collective dynamical behaviors 
have been expected to emerge [|j , . 

On the other hands, the I-V curves along the c-axis in 
Bi-2212 SC's are composed of many branches with almost 
equal inter-branch spacging and the number of those is 
equal to the number of junctions Q, It has been 
widely believed that such structures manifests the inde- 
pendence of junctions inside IJJ's, that is, the coupling 
between junctions is negligible, and then, each junction 
behaves independently [||. Thus, such a understanding 
for the I-V curves is irreconcilable with the above in- 



terpretation for the MRA experiments, and the under- 
standing for IJE's is now still controversial for the two 
experimental results. 

In this paper, we explain the DBCN inside the su- 
perconducting layer in IJJ's, and give a proper model 
Lagrangian for IJJ's. By using the Lagrangian we prove 
that the c-axis propagating longitudinal plasma mode ex- 
ists, and show that the MBS in the I-V characteristics can 
be well reproduced by careful numerical simualtions for 
the model equation derived from the Lagrangian 0, [1). 
Thus, we clarify that both the MRA and the MBS have 
the common origin, i.e., the DBCN. 

Now, let us explain about the DBCN in IJJ's. Fig- 
ure. 1(a) and (b), respectively, are a schematic figure for 
physical events followed by the electron tunneling in con- 
ventional stacked SIS junction (CSISJ) system and IJJ. 
From the figures, it is found that the charges are confined 
near the surface in the CSISJ system while the whole 
of the superconducting layer is charged in the IJJ, i.e., 
the DBCN occurs. This is because the superconducting 
Cu02 layer in IJJ's is extremely thin (~ SA ), i.e., the 
thickness is comparable to the charge screening length 
and therefore the charge screening becomes incomplete 
within the single superconducting layer. In CSISJ sys- 
tems, the electric field generated at a junction site is com- 
pletely screened out by sufficiently thick superconducting 
layers, and the coupling between different junction-sites 
is not formed. On the other hands, in IJJ's, the electric 
fleld generated inside a junction site affects the neigh- 
boring junction sites and therefore the coupling between 
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junctions emerges 0,1^. 
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FIG.l Schematic views of (a) the conventional stacked 
SIS junction (CSISJ) and (b) the intrinsic Josephson 
junction. The arrows indicate the direction of the elec- 
tron tunneling and the generated electric field. In the 
CSISJ, the charges are screened within the range of the 
charge screening length as indicated by the horizontal 
dotted lines, while in IJJ's the whole superconducting 
layer is charged due to the incomplete screening within 
the single superconducting layer. Consequently, the elec- 
tric field generated in a junction site penetrates into the 
neighboring junction sites as the dotted arrows. 

Here, let us give a proper Lagrangian, which describes 
the above DBCN in IJJ's, 



n . ( 



1 — cos Pi 



(1) 



where s and D, respectively, are thickness of the super- 
conducting and insulating layer, /x is the charge screening 
length, e is the dielectric constant, E'^+i,^ is the elec- 
tric field (_L layers) between {i + l)th and ith supercon- 
ducting layer, e* = 2e, and Pi+i^i{= 9g+i{t) - 0e{t) - 

J- J^^^^''^ dzAz{z,t)) is the gauge invariant phase dif- 
ference. In eq.(l) the charge density of £-th supercon- 
ducting layer, pi{t), is related to the scalar potential as 
/^^(*) = ~T^(^o{z£) + -^dtOi) in the gauge-invariant 
form. This Lagrangian corresponds to a discrete ver- 
sion of the well-known time-dependent Ginzburg-Landau 
theory at T = OK for stacked multi-Josephson- 

junction systems when parameters /x and jc are taken as 

H = ^m*vj/{12TTe*'^A'^) and jc = he* A'^ / {M D'^) , where 

Vf is the Fermi velocity, A is the amplitude of the order 
parameter at T = OK, and m* and M are the effective 
masses parallel and perpendicular to the layers, respec- 
tively. The first term in eq.(l) represents the interac- 
tion between the charge density and the gauge invariant 
scalar potential. In conventional theories for CSISJ sys- 
tems, this term is dropped and the dynamics originating 



from only the 2nd and 3rd terms in eq.(l) has so far 
been intensively studied This is because the 

coefficient s//x^ is considered to be very large in CSISJ 
systems and, thus, £ is assumed to have a deep minimum 
at dtOi — —{e* /TijAoizg), which leads to the Josephson 
relation. Thus, the deviations around this minimum do 
not give low energy excitations and the Josephson rela- 
tion always holds within the low energy fiuctuations in 
terms of the phase dynamics in such the CSISJ systems 
Ip^ . However, this is not the case in HTSC's, since the 
thickness, s, of the double superconducting Cu02 planes 
is comparable with ^. This means that the fluctuations 
arising from the first term in eq.(l) cannot be neglected 
in HTSC's, that is, the energy scale of the charge fluctu- 
ations inside the superconducting layers along the c-axis 
is comparable with that of the AC Josephson effect in 
HTSC's. This is the essential point which discriminates 
the IJE's from the AC Josephson effects in conventional 
serial Josephson-junction arrays. 

Let us now prove that the system described by the La- 
grangian (1) has the longitudinal plasma mode. To study 
the response to a longitudinal electric field along the c- 
axis, we derive the dielectric function in the present sys- 
tem. As is well known, in charged systems the Goldstone 
mode is absorbed into the longitudinal gauge field and the 
gauge field becomes massive ( the Anderson-Higgs mech- 
anism ). In order to describe this situation correctly, it 
is convenient to utilize the phason gauge proposed by 
Matsumoto and Umezawa in which the gauge con- 

\i|+Mzil±M =0, where 



dition is given by dtAo{ze) + -f- ^ 
Af^-^^j = JI^^^^^ dzAz{z)/D, and vb is the velocity of 

the phason mode and is given by \J ^"^if {jcY^jcD- 
Here, note that the neutral version of Lagrangian (1), 
in which dynamics of the gauge field Aq[zi) and Afj^-^ ^ 
are dropped, gives the Euler equation for 9( in the linear 

0, where A'^) 
For the gauge 
AUi^t + {hc/e*D){xt+i-Xi). 
Ao{zt) Ao{zi) — {h/e*)dtXii the phason gauge condi- 
tion imposes the equation for xi^ [~'9? + ^b7J2A^^^]x^ = 
0, which is the same as the Euler equation for Op, in the 
neutral case. This result implies that the phason mode 
can be eliminated by the gauge transformation in the 
phason gauge, and the Lagrangian (1) is rewritten in the 
phason gauge as follows. 



approximation as [—9^ -I- u^-^ji A^^^] 6'^ 
is defined as A^^)^^ = Op+i - 29 p 



transformation, A 



'-Jc(l-COS— 



'-^EUi}- (2) 



Variation of for eq. (2) yields one of the Maxwell equa- 
tions, 



El+i l — Egj^i — :^Aq{zi). 



(3) 



2 



From the relation (3), -Bf+i,^ = —dtAg_^_-^^ g/c—[Ao{zi+i) — 
Aq(z£)]/ D, and the phason gauge condition, by ehminat- 
ing the vector potential, we also have another relation 
between the electric field and the scalar potential. 



^A(^))Ao(.,). 



(4) 



Then, from eqs.(3) and (4) one finds the dielectric func- 
tion for the electric field perpendicular to the junctions, 



e{uj, kz) = 1- 



- 2acj2i[l - cosk;,{s + D)] 
In deriving eq.(5) we used the relations. 
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and 



(5) 



As 



eh - exl - ^pi and = ^irJc 
seen in eq.(5), the dielectric function has zero points at 
LLj(kz) — uJpl^y^ + 2a[l — cos kz{s + D)], and the longitu- 
dinal plasma mode propagating along the c-axis is found 
to exist. Here, it should be noted that when the first term 
in the Lagrangian (1) is dropped the obtained Josephson 
plasma mode has no dispersion. 

Next, let us examine dynamics of the Euler-Lagrange 
equation for the gauge invariant phase difference Pi^ij{t) 
derived from the Lagrangian (1) in order to investigate 
the resistive states in IJJ's. The equation is given by two 
Maxwell equations derived from the variation of Aq and 
Az for the Lagrangian and the modified Josephson rela- 
tion obtained by using the gauge invariant charge density 
expression, dtPi+i,i{t) - ^Vi+ij{t) - ^n^{^p^^^(t) - 
Pe{t)), as follows, 



\d'iPi+i,t{t) + sinP^+i,f (t) + —dtPi+i.i{t) = - 

^pl ^pl Jc 

- a (sin Pi+2/+i{t) - 2 sin Pi+i/{t) + sin Pij-i{t)) , (6) 

where the transport current / and the dissipative 
(quasi-particle) current are introduced and those sat- 
isfy the current conservation relation, ^(^jcSvaOi^i^i + 

aEi+i,e^+^dtEi+i^i = /. Also, Upi = and /? = 

^^^^(= where /3c is the McCumber parameter, 
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and a(= ^jj) gives the coupling between junctions. It is 
clearly found that the parameter a can not be neglected 
only in IJJ's due to its definition. This model equation 
has been phenomenologically given by two of us (T.K. 
and M.T. ) Q, and afterwards, Ch.Preis et al. confirmed 
that the same equation can be derived in their leading or- 
der Here, note that if the first term of Llzj is fixed 
to be zero, the model equation becomes the RCSJ model 
[ p2[ , i.e., the model equation for independent junctions. 
Also, it should be noted that when / and p are assumed 
to be zero the linearized eq. (6) has a plane wave solution 
with the same dispersion relation as the relation derived 
from eq.(5) 0, ^. In order to obtain the MBS, we nu- 
merically solve eq.(6) for ten identical junctions under 



the periodic boundary condition. The values of parame- 
ters in eq.(6) a and /3, respectively, are chosen as 0.1 and 
0.20, supposing Bi-2212 case. Although in our previous 
paper ||9[| a was taken as 2.27, the present value employed 
was proved to be reasonable in Bi-2212 ]l^ . 

Figure 2(1) is a hysterisis loop in I-V curve for eq.(6). 
In obtaing this curve, the current I is increased first up to 
a value above the critical current jc and then decreased to 
zero along the arrows in Fig. 2(1). In the current increas- 
ing process, we find a jump at jc, where all junctions syn- 
chronously switch into same resistive states. In Fig. 2(11), 
the time developments of the normalized Josephson cur- 
rents sinPi+i^^ at the point B are plotted for all junction 
sites. As seen in this figure, all junctions show a syn- 
chronous whirling motion. In decreasing the current from 
the point B, we find several step-like structures as seen 
in Fig. 2(1). Figure 2(111) and (IV) show the time devel- 
opments of sinPi^ix of all junction sites at the point C 
and D, respectively. It is found that some junctins are in 
the resistive state, in which sinP^j^u periodically varies 
from -1 to 1, while the others are in the superconduct- 
ing state, in which it shows tiny oscillatory motions. By 
monitoring the time evolution of sinPe+i^e in the region 
around the point C and D as Fig. 2(111) and (IV), one 
may see that in decreasing the current the number of the 
resistive junctions decreases together with the step-like 
structures as seen in Fig. 2(111) and (IV). Such the be- 
haviors are never observed in the case without the cou- 
pling, i.e., a = in eq.(6), which shows only a switching 
behavior into the superconducting state of all resistive 
junctions at a constant current value. 
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FIG. 2 (I) The I-V hysterisis curve obtained by the cur- 
rent increasing and decreasing along the arrows. The 
time developments of sinPi-\-i/ of all junctions bom£ = 
to = 9 at (II) the point B, (III) the point C, and (IV) 
the point D in the I-V hysterisis curve (I). 

Let us now study the region around the point C and 
D in more details. In Fig. 3(a), we present the I-V curve 
enlarged in the region around C and D. In this region, 
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we can get different I-V curves when intervals 51 of the 
decrease in the current value is changed. In the Fig. 3(a), 
two cases with different 61 are plotted. Those behaviors 
imply that several bifurcation points are distributed in 
a very narrow region for the controlling parameter / of 
eq.(6). This is a contrast to the case of a = 0. Thus, 
in order to obtain all possible branches in trajectories of 
eq.(6) in the current decreasing process, we perform the 
calculations many times, changing the decreasing current 
interval SI. Then, when the current is reincreased from 
suitable points of the obtained branches along the arrows 
as seen in Fig. 3 (a), we can obtain several I-V branches, 
which finally lead to the MBS as shown in Fig. 3(b). The 
simulated MBS is composed of equidistant branches and 
the number of resistive branches (=10) is equal to the 
number of junctions (=10). However, in the present sim- 
ulation, it should be noted that we introduce Nyquist 
noise to get the first and ninth branches. But, after the 
transitions by the introduction of the noise, those two 
branches could be stably traced without the noise in the 
current reincreasing. 
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FIG. 3 (a) The enlarged view of the transition region 
in I-V curve. The circles and square points show differ- 
ent TV curves obtained by changing the interval of the 
decrease in the current value, (b) The MBS in the I-V 
curve obtained by the current reincreasing process along 
the arrows as indicated in (a). 

Thus, all the resistive branches can be understood to 
correspond to stable trajectories in the phase space of 
eq.(6), and the MBS as observed in Bi-2212 SC's can be 
reproduced by eq.(6) without inhomogeneities of junc- 
tion parameters. In addition, we mention that if a free 
boundary condition, in which the top and bottom junc- 
tions are distinguished from internal ones, is employed, 
we can obtain all the branches without any noise. In the 
actual experiments for the MBS, the current decreasing 
and increasing are perfomed many times. The present 
simulations also require a number of careful current de- 
creasing and increasing process. 

Finally, let us consider why eq.(6) can give MBS's as 
observed in Bi-2212. The question resolves itself into the 
following two points. The first one is why several step- 



like transitions occur in decreasing the current as seen 
in Fig. 3(a) and the second one is why those branches 
can keep their stablilities in the current reincreasing pro- 
cess as seen in Fig. 3(b). In this paper, we focus on only 
the first point. This is because Takeno et al . have al- 
ready touched the second point in a context for stabil- 
ities of localized whirling (rotating) motions in the dif- 
ferent coupled rotator models with the coupling given 
by trigonometric functions [|l5| and they proved that the 
bounded character of the trigonometric function is es- 
sential for the existence of the localized rotating motions 
JlSf . As shown before, we observe that those step-like 
jumps in the current decreasing result from the step- 
wise transitions of resistive junctions into the supercon- 
ducting state as seen in Fig. 2 (III) and (IV). This im- 
plies that although all junctions are identical each junc- 
tion can switch into the superconducting state at a dif- 
ferent current value I/jc- These behaviors closely re- 
semble those of the independent RCSJ model under the 
presence of inhomogenious distributions in junction pa- 
rameters. We easily notice from the right hand side 
of eq.(6) that the total driving force ( current ) act- 
ing on a junction can be classfied into two components, 
which are the normalized external current I/jc and the 
normalized Josephson current in neighboring junctions, 
Inbj = — Q^(sinPf-|-2.£+i(i)+sinPf^£_i(t)), where the term 
2asmP£+i,^ is transposed to the left hand side. Thus, it 
is clearly found that the total driving force acting on 
each junction depends on the dynamical state of neigh- 
boring junctions. Here, in order to demonstrate it, let us 
suppose that ^-th junction is the resistive state and con- 
sider two simple cases, in which both of the neighboring 
{£± l)-th junctions are in the superconducting state and 
both are in the resistive one. In the former case, since 
the neighboring junctions show only the tiny oscillations 
around the stable points, Inbj is very small ( Inbj ^ ), 
and the total driving current almost equal to the external 
current I/jc- On the other hands, in the latter case, the 
minimum of the instanteneous driving current is given 
by I/jc — 2a, and it can sufficiently become smaller than 
that of the former case. This clearly implies that in the 
current decreasing process the superconducting transi- 
tion of ^-th junction in the latter case occurs prior to the 
former case. Thus, it is found that the switching cur- 
rent value of each junction is dependent on dynamical 
states of neighboring junctions and therefore the super- 
conducting transition points in eq.(6) are distributed in 
the range of the controlling parameter /. From those 
considerations, it is found that the origin of the dynami- 
cal transitions between branches is the coupling between 
junctions due to the charge incomplete screening within 
the single superconducting layer, i.e., the DBCN inside 
the superconducting layer. 

In conclusion, we show that the DBCN inside the su- 
perconducting layer is essential for atomic-scale stacked 
Josephson junction systems as IJJ's, and give the proper 
model Lagragian for IJJ's describing the DBCN. We clar- 
ify that the Lagragian involves the dynamics of the lon- 
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gitudinal propagating Joscphson plasma responsible for 
the MRA observed in the longitudinal configuration while 
the Euler-Lagrange equation under the transport and the 
disspative current can reproduce the MBS. 
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